Birutín

Introducción

Aunque sean estructuras sencillas, la vida de un virus en absoluto lo es. No puede ser demasiado letal ni demasiado inocente, tiene que aprender a vivir con su hospedador sin exprimirlo mucho ni poco, en un juego evolutivo con múltiples facetas e implicaciones..

Lo que sigue a continuación no es más que un simple ejercicio en el que se pone a prueba el poder de predicción de unas simples ecuaciones y un juego con sencillas reglas. Sus implicaciones, por tanto, son más heurísticas que aplicables. Cualquier modelo que verdaderamente pueda ser predictivo para el caso que todos tenemos en mente debe considerar muchos aspectos secundarios que aquí se obvian.

No obstante, los resultados son lo bastante interesantes como para hacernos pensar sobre detalles aparentemente triviales que después se manifiestan de forma sorprendente.

Matemática de una epidemia

Una epidemia es el fenómeno demográfico provocado por la dispersión de una enfermedad contagiosa en una población. La enfermedad va saltando de individuo en individuo a través de un vector, un virus en este caso. La epidemia es la dinámica temporal que manifiesta el número de infectados, que suele comenzar de una manera explosiva (exponencial, propiamente hablando) para llegar a un pico y luego disminuir hasta que desaparece la enfermedad de la población.

Ecuaciones

El modelo más sencillo surge a partir de los siguientes axiomas:

• Cualquier individuo susceptible puede contraer la enfermedad con cierta probabilidad en un plazo determinado.

• Esta probabilidad aumenta proporcionalmente al número de infectados.

• Al cabo de cierto tiempo la enfermedad cesa en el infectado, bien con su inmunización o con su muerte, con el efecto neto de sacarlo de la población susceptible o infectada.

Con estas premisas de partida, el par de sesudos matemáticos ingleses Kermack y Kendrick a partir de unos trabajos previos muy enrevesados del teniente coronel Ross y (probablemente) su querida madame Hudson, desarrollaron hace un siglo un completa teoría matemática que, simplificada al considerar tasas constantes y una población suficientemente grande, se puede resumir en la siguiente ecuación:

\[\begin{align} cambio\;en\;infectados =& \;contagios - recuperaciones \\[7pt] \frac{dI}{dt} =& \beta S I - \gamma I \end{align}\]

El término \(\beta S I\) representa el número de nuevos contagios por unidad de tiempo, siendo \(S\) e \(I\), respectivamente, el número de susceptibles e infectados en la población; mientras que \(\beta\) es un coeficiente de contagio que aúna en un solo parámetro las propiedades contagiosas de la enfermedad y los hábitos sociales de la población: cuanto mayor es, más fácilmente se contagia la enfermedad.

Por otro lado, \(\gamma I\) representa la retirada de individuos infectados (bien porque se recuperan volviéndose inmunes, bien porque pasan a mejor vida; el caso es que ya no juegan). El coeficiente \(\gamma\) está relacionado con la duración de la enfermedad, en particular, con el tiempo que el individuo es contagioso, como veremos con detalle dentro de un rato.

Si incluimos la dinámica de los susceptibles \(S\) y los retirados \(R\), tendremos esta terna de ecuaciones:

\[\begin{align} \frac{dS}{dt} &= - \beta S I \\[7pt] \frac{dI}{dt} &= \beta S I - \gamma I \\[7pt] \frac{dR}{dt} &= \gamma I \end{align}\]

Ejemplo académico

Para comprender mejor el sentido de ecuaciones y parámetros, veamos cómo aparece el modelo a partir de unas observaciones naturales y más intuitivas. Imaginemos que, en cierta población con 50.000 habitantes susceptibles ocurre una epidemia y se ha comprobado que:

• Un infectado contagia la enfermedad a 3 personas, en promedio.
• La enfermedad en fase contagiosa dura 1 semana, en promedio.

Cálculo de \(\beta\)

Vamos con la primera observación. ¿Podemos deducir el valor de \(\beta\)? A ver…

1. Un infectado contagia a 3 susceptibles; como lo puede hacer mientras dura su fase contagiosa (7 días), cada día contagiará a \(\frac{3}{7}\) individuos.
2. Como hay \(I\) infectados, el número de nuevos contagios diarios será \(\frac{3}{7}I\).
3. ¿Y cómo metemos \(S\)? Veamos… La primera observación es claramente incompleta, le falta «…al principio de la epidemia»: está claro que la probabilidad de que un infectado se lo pegue a un susceptible es mucho mayor al principio de la epidemia (cuando \(S\) es grande) que al final (cuando \(S\) es menor), por pura ley de acción de masas. Esto significa que \(\frac{3}{7}\) es justamente \(\beta S_0\), siendo \(S_0\) la población susceptible inicial (50.000). Ya lo tenemos:

\[\beta = \frac{3}{7}·\frac{1}{50000} = 8.57 \cdot 10^{-6} \; (\text{individuos}^{-1}\text{·día}^{-1})\]

Cálculo de \(\gamma\)

Vamos con la segunda. ¿Podemos deducir el valor de \(\gamma\)? Para verlo más fácilmente, supongamos que tenemos 35 infectados; ¿cuántos se recuperarán en un día? La figura tal vez ayude:

Figura 1. Significado de \(\gamma\). Como la duración de la enfermedad es de una semana, cuando pase 1 día se recuperarán 5 de los 35 infectados, es decir \({1}/{7}\) que es justo el valor de \(\gamma\). Así que \(\gamma\) es, intuitivamente, el inverso de la duración de la enfermedad.

En principio, porque ya veremos más adelante. El razonamiento anterior asume que tenemos infectados recientes y veteranos en la misma proporción (los puntitos rojos se distribuyen homogéneamente en el plazo de convalecencia). No es difícil demostrar que esta premisa se cumple siempre en un estado estacionario pero, como la epidemia cursa en forma de onda, esto supondrá un problema que solucionaremos más adelante.

Dinámica

Ya tenemos todo. Si integramos las ecuaciones veremos el aspecto gráfico de la epidemia:

{
# INTEGRACIÓN DE LAS ECUACIONES S.I.R.

#---- tiempo final
tiempo_final <- 100
t <- seq(from = 1, to = tiempo_final)

#---- condiciones iniciales y valores de coeficientes
S <- numeric()
S[1] <- 50000 # susceptibles iniciales
I <- numeric()
I[1] <- 1      # infectados iniciales
R <- numeric()
R[1] <- 0      # retirados iniciales
beta <- 8.57E-6
gama <- 1/7

#---- bucle integración Euler
for (i in 1:(length(t)-1)){
  contagios <- beta*S[i]*I[i]
  retirados <- gama*I[i]
  S[i+1] <- S[i] - contagios
  I[i+1] <- I[i] + (contagios - retirados)
  R[i+1] <- R[i] + retirados
}

#---- gráfico
plot(t, S, ylim = c(0, 50000), type = "l", col = "green", lwd = 2,
     xlab = "TIEMPO (días)", ylab = "INDIVIDUOS")
lines(t, R, col = "orange", lwd = 2)
lines(t, I, col = "red", lwd = 2)
legend(0, 50000, c("S","I","R"), lwd = 2, col = c("green", "red", "orange"), cex = 0.8, bty = "n")
}

Figura 2. Dinámica de una epidemia. La línea roja representa la epidemia, propiamente dicha. Puede observarse el inicio exponencial, el pico y la ulterior decaída, fases típicas en toda epidemia. Ojo a la línea verde: no todos los susceptibles deben pasar la enfermedad antes de que cese la epidemia.

Características de la dinámica epidémica

En la Figura 2 podemos apreciar las características dinámicas de toda epidemia: el inicio exponencial, el pico de infectados y el descenso final hasta la desaparición de la enfermedad en la población. Un aspecto interesante es que la epidemia puede (aunque no siempre) acabar antes de que todos los susceptibles pasen la enfermedad.

Otra característica importante es la asimetría de la curva: el inicio es mucho más empinado y el descenso más lento. Esto ocurre porque el contagio depende del cuadrado de los individuos mientras que la recuperación es función lineal.

Efecto de la densidad de población

El tamaño de la población también es un factor crítico. Podemos verlo mejor si expresamos los infectados en % de la población y hacemos varias simulaciones con distintas densidades de población, por ejemplo con 400.000, 200.000, 100.000, 25.000 y 12.000 almas:

Figura 3. Efecto de la densidad de población. La misma enfermedad desencadena distintos tipos de epidemia dependiendo del tamaño de la población susceptible: cuanto mayor es la población, más rápidamente se desencadena la epidemia, mayor porcentaje de la población se afecta y antes acaba el proceso. Pero… ¿qué ocurre con la de 12.000 habitantes?

Efectivamente, cuanto mayor es la población, más devastadora es la epidemia. Los efectos en poblaciones menores son cada vez más suaves, llegando al extremo, en la ciudad con 12.000 habitantes, de ni siquiera producirse la epidemia, a pesar de tratarse de la misma enfermedad y el mismo tipo costumbres sociales. ¿Por qué?

La mayor virulencia en poblaciones grandes podemos comprenderla intuitivamente: si un infectado contagiaba a 3 personas en una ciudad de 50.000 individuos, entonces, en una ciudad de 400.000 habitantes (8 veces más poblada) contagiará a 24 (8 veces más) y, por tanto, la velocidad de propagación será 8 veces mayor que en la ciudad de 50.000 almas.

El factor \(R_0\)

Al factor por el que se multiplica la infección (o número de individuos a los que contagia un infectado) se denomina \(R_0\) y vemos que es directamente proporcional a la densidad de población. Podemos calcular su valor sin necesidad de acudir a las diferenciales: al invertir el razonamiento mediante el que calculamos \(\beta\) y \(\gamma\) vemos que

\[R_0 = \frac{\beta S}{\gamma}\]

Ya tenemos su valor para la ciudad de 50.000 personas (era uno de los axiomas) y hemos hecho el cálculo para la población de 400.000 habitantes. Por curiosidad, apliquémoslo a las otras poblaciones:

HABITANTES	\(R_0\) (al inicio de la epidemia)
400.000	24
200.000	12
100.000	6
50.000	3
25.000	1.5
12.000	0.75

¡Ah migo! En la última línea tenemos la explicación del fracaso de la epidemia en la ciudad pequeña: si un infectado solo es capaz de transmitir la enfermedad a 3/4 de persona, está claro que la epidemia no se producirá. Para que explote la epidemia, \(R_0\) debe ser mayor que 1.

Como depende directamente de \(S\) y éstos van siempre disminuyendo, el factor \(R_0\) también va disminuyendo desde el inicio de la epidemia. Veámoslo en el primer gráfico:

Figura 4. Evolución de \(R_0\). El mismo gráfico de antes (50.000 susceptibles), pero ahora con el valor de \(R_0\) superpuesto. Su curva es idéntica a la de S, de hecho, es la misma multiplicada por \({\beta}/{\gamma}\). Observa que, justo cuando \(R_0\) cruza el umbral de 1, la epidemia alcanza su pico y comienza a declinar.

Ejemplo real

El inicio de la epidemia de COVID en la Costa del Sol nos da un ejemplo, tan ilustrativo como triste, de estas simple matemáticas. Las primeras semanas de propagación de la enfermedad ocurrieron con esta cifras oficiales:

FECHA	INFECTADOS
27/2/20	1
28/2/20	4
29/2/20	6
1/3/20	6
2/3/20	6
3/3/20	7
4/3/20	7
5/3/20	9
6/3/20	18
7/3/20	22
8/3/20	27
9/3/20	42
10/3/20	59
11/3/20	77
12/3/20	97
13/3/20	128
14/3/20	180

Figura 5. Primeras semanas de la COVID en la Costa del Sol. El aumento de infectados se produjo de forma exponencial.

Cálculo de los parámetros en Málaga

Veamos si podemos calcular los parámetros de la epidemia. El problema es muy interesante, matemáticamente hablando, pues consiste en el ajuste de una función que no tiene una expresión analítica explícita, pues resulta de la integración numérica de dos ecuaciones diferenciales. Dejo un frente abierto, intentando ajustar tanto \(\beta\) y \(\gamma\) como \(S_0\) de forma independiente, porque ahora el objetivo principal es más evolutivo que matemático. No obstante, con un poco de habilidad se puede simplificar el asunto, si podemos estimar el valor de \(S_0\) y \(\gamma\) independientemente de los datos.

Valor de \(S_0\)

Propiamente dicho, no es un coeficiente, pero sí un parámetro del modelo que debemos ajustar. Dado que el foco principal se localizó en toda la Costa del Sol (1.400.000 habitantes) y esta zona muestra propiedades de movilidad y costumbres muy homogéneas, parece acertado empezar con este valor.

Valor de \(\gamma\)

La duración total de la enfermedad es muy indefinida, va desde los 24 a los 57 días, según las fuentes de donde proceda la información. No obstante, dado que en Málaga ya se aislaban a las personas con los primeros síntomas, para nuestro cálculo solo tendremos en cuenta el tiempo que el recién contagiado circuló libremente sin mostrar síntomas pero siendo capaz de contagiar a otros, un periodo que puede ser de unas dos semanas. Tomaremos un valor de 15 días y luego calcularemos la incertidumbre del resultado con respecto a este dato, si es menester. Así, tomaremos un valor de \(\lambda = 0.07 \text{ días}^{-1}\).

Valor de \(\beta\)

La ecuación central del modelo se puede escribir de esta forma:

\[\frac{dI}{dt} = (\beta S - \gamma) I\]

como aún aparece \(S\) en la expresión, su ajuste sigue requiriendo el método más sofisticado que tal vez abordemos en otra ocasión. No obstante, en aquel tiempo el número de infectados era muy pequeño (en comparación con \(S_0\)) y podemos asumir que \(S\) es prácticamente constante en estos primeros compases, de manera que:

\[\frac{dI}{dt} \approx (\beta S_0 - \gamma) I\]

expresión que tiene integral inmediata y da la siguiente ecuación:

\[I = I_0 e^{(\beta S_0 - \gamma) t}\]

cuyo ajuste es inmediato por linealización. Sin embargo, para evitar el sesgo que esto provoca sobre los números más grandes, vamos a ajustar la expresión directamente por optimización:

{
#---- datos reales
I_real <- c(1,4,6,6,6,7,7,9,18,22,27,42,59,77,97,128,180)
n <- length(I_real)
t_real <- seq(1:n)
plot(t_real, I_real, xlab = "TIEMPO (días)", ylab = "INFECTADOS")

#---- parámetros estimados
S0 <- 1400000
gama <- 0.07

#---- definir la suma de diferencias cuadraticas como función
SDC <- function(dat, par) {
  I <- par[1]*exp((par[2]*S0 - gama)*t_real)
  with(dat, sum((I_real - I)^2))
}

#---- buscar los parámetros que la minimizan
datos <- data.frame(t_real,I_real)
resultado <- optim(par = c(1, 1E-7), SDC, dat = datos)
I0 <- resultado$par[1]
beta <- resultado$par[2]

#---- gráfico de la nueva curva ajustada
t <- seq(0, 20, 0.1)
I <- I0*exp((beta*S0 - gama)*t)
lines(t,I, type = "l", col = "red", lwd = 2)

#---- cálculo del coeficiente de determinación de la nueva curva
I_aj <- I0*exp((beta*S0 - gama)*t_real)
R2curva <- 1 - sum((I_real - I_aj)^2)/((n-1)*var(I_real))
text(12,150,bquote(R^2:.(R2curva)))

#---- información sobre beta e I0
text(5,150,bquote(beta:.(beta)))
text(5,130,bquote(I_0:.(I0)))
}

Figura 6. Ajuste de \(\beta\). A pesar de haber estimado dos parámetros y ajustado el tercero, este simple modelo es capaz de explicar más del 99% de la varianza de los datos (\(R^2 = 0.997\)) y nos da un valor de \(\beta =2.64·10^{-7}\) individuos⁻¹ día⁻¹.

En este modelo ajustado, tenemos un valor de

\[R_0 = \frac{\beta S_0}{\gamma} = \frac{2.64·10^{-7} · 1400000}{0.07} = 5.28\]

que coincide con los cálculos hechos por modelos más sofisticados y está en concordancia con los valores de \(R_0\) determinados inicialmente en China y, durante bastante más tiempo, en Italia.

Como curiosidad

De esta sencilla aproximación se puede deducir una conclusión de interés epidemiológico: para reducir el \(R_0\) de la COVID por debajo de \(1\), dado que es muy difícil acelerar \(\gamma\), se puede disminuir \(S\) (con vacunas) o reducir \(\beta\) (con restricciones de movilidad). Según el valor calculado arriba, bastaría con dividir algunos de estos dos elementos entre \(5.3\), esto es, en un 82 % aproximadamente. Esto significa vacunar al 82 % de la población susceptible, o bien, limitar nuestros contactos a 1 de cada 6 personas con las que tratamos habitualmente. Curiosamente, con solo los datos de las primeras semanas y un sencillo análisis, hemos alcanzado las mismas conclusiones que ahora guían la lucha contra esta epidemia… ¡Seguro que ha sido casualidad!

Epidemia con reinfección

Cuando los recuperados pierden la inmunidad y se vuelven suceptibles al cabo de cierto tiempo, aparece una dinámica cíclica conocida como olas. Esto ocurre porque se introduce un nuevo elemento en las ecuaciones que las convierte en un oscilador amortiguado:

\[\begin{align} \frac{dS}{dt} &= \alpha R - \beta S I \\[7pt] \frac{dI}{dt} &= \beta S I - \gamma I \\[7pt] \frac{dR}{dt} &= \gamma I - \alpha R \end{align}\]

El nuevo coeficiente \(\alpha\) es inversamente proporcional al tiempo que dura la inmunidad de un recuperado (ya no pueden ser muertos, claro está).

Veamos la dinámica del primer ejemplo, suponiendo que la inmunidad dura 20 semanas (\(\alpha = 1/140\) días⁻¹)

Figura 7. Epidemia con reinfección. Cuando los recuperados pierden la inmunidad se producen oleadas de infectados hasta que cesa la epidemia. La distancia entre los picos de estas olas y su intensidad son proporcionales a la duración de la inmunidad. Los coeficientes son los mismos de la Figura 2, más una pérdida de inmunidad a los 140 días.

También se podría añadir una tasa de mortalidad asociada con la enfermedad, para que el asunto fuese más realista, aunque el interés matemático es mínimo pues solo se manifiesta en una caída exponencial de los recuperados y la dinámica oscilante seguiría siendo prácticamente la misma.

Modelo basado en individuos

Lo de antes era solo para recordar lo básico del asunto. Aquí empieza lo interesante.

Cuando se afronta el fenómeno epidémico desde la mecánica basada en individuos se disfruta de unos axiomas más intuitivos:

• Los individuos se mueven al azar por su pueblo.

• Cuando un infectado contacta con un susceptible lo puede contagiar con una probabilidad \(p\).

• Al cabo de cierto tiempo, el infectado se recupera y ya no contagia más.

  Aquí tenemos propiedades individuales, en vez de las colectivas del modelo diferencial. No nos enrollemos más y veamos la película del movimiento individual y la transmisión de la enfermedad:

Vídeo 1. Dinámica espacial de la epidemia. Los puntos rojos representan a los infectados (como se esperaba) que, cuando entran en contacto con los susceptibles (verdes) los infectan con cierta probabilidad. Al cabo de algunos días, los infectados se recuperan (negros).

Figura 8. Simulación basada en individuos. Con los titubeos debidos al azar, el patrón de la epidemia aparece claramente.

Aquí debajo el código:

{
# MODELO S.I.R. BASADO EN INDIVIDUOS
# Modelo básico para explicar el desarrollo de una epidemia
# (sin vídeo)


#---- tiempo final y tamaño del pueblo
tiempo_final <- 150
xmax <- 60                    # ancho del pueblo (desde x = 1 hasta x = 100)
ymax <- 60                    # largo "     "    (desde y = 1 hasta y = 100)
XY <- matrix(0L, xmax, ymax)   # el pueblo, propiamente dicho

#---- variables
St <- integer()
It <- integer()      # para guardar los valores de S, I, R a cada tiempo
Rt <- integer()

x <- integer()      # posición x de los individuos
y <- integer()      #     "    y      "       "
estado <- integer() # estado del individuo  S:0  I:1  R:2
conval <- integer() # tiempo que lleva de convalecencia

#---- coeficientes y demás
p <- 0.9                 # probabilidad de contagio
d <- 15                  # duración de la enfermedad
mov <- c(-1, 0, 1)       # movimiento (una casilla pacalao)
pmv <- c(0.4, 0.2, 0.4)  # probabilidad de moverse en cada dirección XY

#---- condiciones iniciales
N <- 1000      # número total de individuos
S <- N - 1     # susceptibles iniciales
I <- 1         # infectados      "
R <- 0         # retirados       "
St[1] <- S
It[1] <- I
Rt[1] <- R

# Localizar a los individuos
for(i in 1:N){
  x[i] <- floor(xmax*runif(1) + 1)  # pensando en C...
  y[i] <- floor(ymax*runif(1) + 1)
  estado[i] <- 0
  XY[x[i], y[i]] <- i              # localizar al último en llegar
}

# infectados
for(i in 1:I){
  x[i] <- floor(xmax*0.5 + 1)  
  y[i] <- floor(ymax*0.5 + 1)
  estado[i] <- 1
  XY[x[i], y[i]] <- i              
  estado[i] <- 1
  conval[i] <- 0
}

#---- BUCLE DEL TIEMPO
for(t in 2:tiempo_final){
  
  XY1 <- XY     # para hacer las operaciones sincronizadas
  XY[] <- 0L    # borrar la matriz
  
  #---- bucle de los individuos
  for(i in 1:N){
    
    # comprobar convalecencia
    if(estado[i] == 1) {# solo para infectados
      conval[i] <- conval[i] + 1
      if(conval[i] > d) {# se recupera
        estado[i] <- 2
        I <- I - 1
        R <- R + 1
      }
    }
    
    # movimiento al azar
    dx <- sample(mov, size = 1, prob = pmv)
    dy <- sample(mov, size = 1, prob = pmv)
    x[i] <- x[i] + dx
    y[i] <- y[i] + dy
    if(x[i] > xmax) {x[i] <- 1} # por si se sale por los lados
    if(x[i] < 1) {x[i] <- xmax}
    if(y[i] > ymax) {y[i] <- 1} # por si se sale por arriba/abajo
    if(y[i] < 1) {y[i] <- ymax}
      
    # detectar encuentro
    o <- XY1[x[i], y[i]]
    if(o > 0) {
      if(o != i) { # se ha encontrado con otro
        # se lo pego
        if((estado[i] == 1) & (estado[o] == 0)){
          if(runif(1) < p){
            estado[o] <- 1
            conval[o] <- 0
            I <- I + 1
            S <- S - 1
          }
        }
        # me lo pega
        if((estado[o] == 1) & (estado[i] == 0)){
          if(runif(1) < p){
            estado[i] <- 1
            conval[i] <- 0
            I <- I + 1
            S <- S - 1
          }
        }
      }
    }
    
    XY[x[i], y[i]] <- i   # actualizar localización
  }
  
  St[t] <- S
  It[t] <- I    # guardar valores
  Rt[t] <- R
  
}

#---- gráfico final
plot(St, type = "l", lwd = 2, col = "green", ylim = c(0,N), xlab = "Tiempo", ylab = "Individuos")
lines(It, lwd = 2, col = "red")
lines(Rt, lwd = 2)
legend(0, N/2, c("S", "I", "R"), lty = c(1, 1, 1), lwd = c(2, 2, 2), col = c(3, 2, 1), bty = "n")
}

No me convence mucho el asunto de las matrices, porque el último en llegar al sitio es el que marca las propiedades. Vamos a dejarlo así por el momento, posiblemente no tenga mucha importancia en pueblos enrarecidos, pero la mosca está en la oreja.

Veamos el promedio de 100 simulaciones

Figura 9. Promedio de 100 simulaciones. Se parece al modelo diferencial peeerooo…

El inicio de la epidemia es demasiado lineal. Podría tratarse de:

1. Explicación sesuda: Aunque haya muchos susceptibles, cada infectado solo puede acceder a los que están más cerca de él, que es una cantidad constante. Esto se traduce en un frente de llama que avanza sobre un mar se susceptibles dejando tras sí un bosque de recuperados.
2. Artefacto: Al usar la misma matriz para todos los individuos, pueden enmascararse infectados por solapamiento con recuperados. Esto se podría evitar con matrices independientes para \(S\) e \(I\) y, ya puestos, también para \(R\).

Antes de ponerme con las matrices independientes, confirmemos que se trata del efecto frente de llama. Para evitarlo, creo que basta con disminuir la probabilidad de contagio y aumentar la duración de la enfermedad. Veamos:

Figura 10. Evitar el frente de llama. Efectivamente, cuando se reduce la probabilidad de contagio y aumenta el tiempo de convalecencia desaparece el efecto frente de llama. Ahora lo sospechoso es la caída simétrica, que no concuerda con lo esperado según un proceso exponencial.

¿Será un efecto tapadera de los recuperados? Si hay muchos, pueden encubrir a los susceptibles. Veamos la simulación quitándolos:

Figura 11. Sin recuperados. Escasa diferencia, así que no merece la pena, porque después el código se volverá más complicado y mucho menos elegante.

¿Será un efecto de la convalecencia, que no tiene un valor constante? Veamos cómo se distribuye la convalecencia:

Vídeo 2. Distribución de tiempo de convalecencia. Conforme van entrando individuos en el periodo de convalecencia, van saliendo por el otro lado. El tiempo medio de convalecencia cumplida (línea roja) no es constante: al principio es menor y, conforme avanza el tiempo, aumenta. Esto significa que no se puede considerar constante y, por tanto, tampoco \(\gamma\).

Figura 12. Recuperación de infectados. El modelo diferencial dispersa en el tiempo la aparición de recuperados al considerar homogéneo el grupo de infectados (como se explicó intuitivamente en la Figura 1). En el Vídeo 2 se ve como la salida de convalecencia se produce según una cola tipo FIFO «primero en entrar, primero en salir». Como consecuencia, la ola de infectados se alarga en el modelo diferencial

Sorprendentemente, al resolver este asunto hemos comprobado que el clásico modelo SIR es realista en la propagación de la epidemia y, sin embargo, falla al explicar el final de la ola porque \(\gamma\), en realidad, no es constante. Por cierto, ¿qué tipo de función sería? ¡No acaban los flecos!